In [2]:
import numpy as np
from matplotlib import rcParams
import matplotlib.pyplot as pl
import scanpy.api as sc
import pandas as pd
import matplotlib.pyplot as plt
import episcanpy.api as epi
import anndata as ad
import numpy as np
sc.settings.verbosity = 3 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.settings.set_figure_params(dpi=100, color_map='viridis') # low dpi (dots per inch) yields small inline figures
sc.logging.print_versions()
#Set working directory
wd = '/path/to/file/sample_pooled_preprocess_revision1/'
direc = wd
sc.settings.figdir = './plots/'
results_file="./write/20190416_snATACseq_embryo_revision01_doublets_cisTopic_50_100_filtered_final.h5ad"
Read files¶
In [3]:
##====== Read counts ======##
filename_data = direc + '11_matrix_afterClusterQC/embryo_revision1_allPeaks_afterClusterPeak.mat.bin.mtx'
filename_gene_names = direc + '11_matrix_afterClusterQC/embryo_revision1_allPeaks_passedQC_peakNames.txt'
filename_barcodes = direc + '11_matrix_afterClusterQC/embryo_revision1_allPeaks_passedQC_barcodeNames.xgi'
#filename_clustersAll = wd + 'data/cellTypes_20180215.txt'
#adata_all_cells = epi.tl.read_ATAC(filename_data, filename_barcodes, filename_gene_names,path_file='')
print('reading counts')
adata_all_cells = sc.read(filename_data, cache=True).transpose()
adata_all_cells.X = adata_all_cells.X.astype(np.int64)
print('reading genes')
adata_all_cells.var_names = np.genfromtxt(filename_gene_names, dtype='str')
print('reading cells')
adata_all_cells.obs_names = np.genfromtxt(filename_barcodes, dtype='str')
In [4]:
adata_all_cells.shape
Out[4]:
colour Palettes¶
In [5]:
#adata = sc.read(results_file)
In [6]:
##====== Create colour palette for gene expression profiles ======##
from matplotlib.colors import LinearSegmentedColormap
rmap = LinearSegmentedColormap.from_list(name='gene_cmap',
colors=['lightgrey', 'thistle', 'red', 'darkred'])
cmap = LinearSegmentedColormap.from_list(name='gene_cmap',
colors=["#BFBFBF","#6495ED","#000000"])
cisTopic with reads_in_peaks > 0.24¶
In [7]:
cells = np.genfromtxt(direc + "12_barcodeStats_celltypePeaks/embryo_revision1_readsPeaks24.xgi", dtype='str')
In [8]:
adata = adata_all_cells[np.array([str(i) in cells for i in adata_all_cells.obs_names]),:]
In [9]:
adata.shape
Out[9]:
In [10]:
from scipy.io import mmread
cistopic = mmread(direc + '13_cisTopic/data/cisTopic_matrix_50_100_afterReadsInPeaksQC_24.mtx')
In [11]:
cistopic = cistopic.todense()
cistopic= cistopic.transpose()
cistopic.shape
Out[11]:
In [12]:
adata.obsm['X_pca'] = cistopic
In [13]:
adata.obs=pd.DataFrame(cistopic)
UMAP visualisation and clustering¶
In [14]:
sc.pp.neighbors(adata)
In [15]:
sc.tl.umap(adata,random_state=4)
In [16]:
sc.tl.louvain(adata,random_state=7)
In [52]:
sc.pl.umap(adata, color=[ 'louvain'])
Plot cistopic scores¶
In [53]:
for i in adata.obs:
sc.pl.umap(adata, color=[i])
Write files¶
In [111]:
#save UMAP
UMAP = pd.DataFrame(data=adata.obsm['X_umap'])
UMAP.to_csv(direc + '14_visualisation/data/snATACseq_embryo_revision1_reads_in_peaks_above_24_UMAP.csv', sep=',')
#save clusters
louvain = pd.DataFrame(data=adata.obs['louvain'])
louvain.to_csv(direc + '14_visualisation/data/snATACseq_embryo_revision1_reads_in_peaks_above_24_Louvain.csv', sep=',')
Count number of cells with regions accessible per cell type¶
In [65]:
mat = np.zeros((len(adata.var_names), len(np.unique(adata.obs['louvain']))))
count=0
for i in np.unique(adata.obs['louvain']):
print(i)
adata_endo = adata[[j == i for j in adata.obs["louvain"]],:]
numEndo = np.sum(adata_endo.X,axis=0).tolist()
print(numEndo[0][0:5])
mat[:,count] = numEndo[0]
count+=1
print(count)
In [66]:
MAT = pd.DataFrame(data=mat)
MAT
Out[66]:
In [115]:
MAT = pd.DataFrame(data=mat)
MAT.to_csv(direc + '14_visualisation/data/snATACseq_embryo_revision01_pooledNumber_OCR_perCelltype.csv', sep=',')
#save clusters
#save peaks used
CELL = pd.DataFrame(data=np.unique(adata.obs['louvain']))
CELL.to_csv(direc + '14_visualisation/data/snATACseq_embryo_revision01_pooledNumber_OCR_perCelltype_celltypeLabel.csv', sep=',')
PEAK = pd.DataFrame(data=adata.var_names)
PEAK.to_csv(direc + '14_visualisation/data/snATACseq_embryo_revision01_pooledNumber_OCR_perCelltype_peakLabel.csv', sep=',')
In [116]:
MAT.shape
Out[116]: